DOES POPULATION DENSITY AND LAND USE AFFECT THE RATE OF REDUCTION OF FUKUSHIMA DAICHI RADIATIONS?

Part I: Data Loading, Cleaning and Exploring

1.1 Nihomatsu City,70km to Fukushima Daichi

1.11 Required Packages

library(leaflet)
library(readr)
library(dplyr)
library(RColorBrewer)
library(Hmisc)
library(ggplot2)
library(formatR)
1.11 Loading June 2011 and Nov 2013 Fukushima Data and selecting Nihomatsu’s.
fuk2011 <- read.csv("jun_2011_fukushima.csv")
fuk2013 <- read.csv("nov_2013_fukushima.csv")
1.12 Change to machine readeable column names
names(fuk2011) <- c("gridcode", "pref", "city", "gridCenterNorthlat", 
    "gridCenterEastlng", "gridCenterNorthlatDec", "gridCenterEastlngDec", 
    "daichi_distance", "no_samples", "AvgAirDoseRate", "NE_nLat", 
    "NE_eLong", "NW_nLat", "NW_eLong", "SW_nLat", "SW_eLong", 
    "SE_nLat", "SE_eLong")
names(fuk2013) <- c("gridcode", "pref", "city", "gridCenterNorthlat", 
    "gridCenterEastlng", "gridCenterNorthlatDec", "gridCenterEastlngDec", 
    "daichi_distance", "no_samples", "AvgAirDoseRate", "NE_nLat", 
    "NE_eLong", "NW_nLat", "NW_eLong", "SW_nLat", "SW_eLong", 
    "SE_nLat", "SE_eLong")
dim(fuk2013)
## [1] 119522     18
dim(fuk2011)
## [1] 45273    18
max(fuk2013$AvgAirDoseRate)
## [1] 39
min(fuk2013$AvgAirDoseRate)
## [1] 0.009
max(fuk2011$AvgAirDoseRate)
## [1] 37
min(fuk2011$AvgAirDoseRate)
## [1] 0.008
1.14 Create air dose quantiles that are plot-able,i.e 6 categorical variables.
fuk2011_q <- fuk2011 %>% mutate(dose_quants = cut2(fuk2011$AvgAirDoseRate, 
    cuts = c(0.1, 1, 5, 15, 25, 35, 45), levels.mean = TRUE))
fuk2013_q <- fuk2013 %>% mutate(dose_quants = cut2(fuk2013$AvgAirDoseRate, 
    cuts = c(0.1, 1, 5, 15, 25, 35, 45), levels.mean = TRUE))

# write.csv(fuk2011_q, file='fuk2011dose.csv',row.names =
# FALSE) write.csv(fuk2013_q,
# file='fuk2013dose.csv',row.names = FALSE)

fuk2013dose <- read.csv("fuk2013dose.csv")
# remove the slash on grides
fuk2013dose$gridcode <- gsub("_", "", fuk2013dose$gridcode)
# remove the last character on grides, 250m to 500m
fuk2013dose$gridcode <- gsub(".{1}$", "", fuk2013dose$gridcode)
  • Visible reduction of Average Air Dose Distribution by half in Nihomatsu.
  • Trouble is knowing the distribution of causes of this reduction?
1.15 Color function
iro <- colorFactor(palette = "YlOrRd", domain = fuk2011_q$dose_quants)
iro2013 <- colorFactor(palette = "YlOrRd", domain = fuk2013_q$dose_quants)
# Link of Daichi
fukulink <- paste(sep = "<br/>", "<br><a href='http://www.tepco.co.jp
                  /en/decommision/index-e.html'>Fukushima Daichi</a></b>", 
    "Source of radiations")
1.16 Fukushima Average Air Dose Rate for 2011
fuk2011_plot <- leaflet() %>% addTiles() %>% addRectangles(data = fuk2011_q, 
    lng1 = ~SW_eLong, lat1 = ~SW_nLat, lng2 = ~NE_eLong, lat2 = ~NE_nLat, 
    color = ~iro(fuk2011_q$dose_quants)) %>% addLegend("bottomright", 
    pal = iro, values = fuk2011_q$dose_quants, title = "AvgAirDoseRate", 
    labFormat = labelFormat(prefix = "µSv/h "), opacity = 1) %>% 
    addPopups(lat = 37.4211, lng = 141.0328, popup = fukulink, 
        options = popupOptions(closeButton = TRUE))
fuk2011_plot


1.17 Fukushima Average Air Dose Rate for 2013
fuk2013_plot <- leaflet() %>%
        addTiles()%>%
        addRectangles(data = fuk2013_q,lng1 = ~SW_eLong, lat1 = ~SW_nLat,
                      lng2 = ~NE_eLong, lat2 = ~NE_nLat,
                      color = ~iro2013(fuk2013_q$dose_quants))%>%
        addLegend("bottomright", pal = iro2013, values = fuk2013_q$dose_quants,
                  title = "AvgAirDoseRate",
                  labFormat = labelFormat(prefix = "µSv/h "),
                  opacity = 1)%>%
        addPopups(lat = 37.4211, lng = 141.0328,popup = fukulink,
                  options = popupOptions(closeButton = TRUE)) 
fuk2013_plot


1.18 Nihommatsu 2011, Counts of Measuring Locations per Air Dose Rate
ggplot(fuk2011_q, aes(daichi_distance,AvgAirDoseRate)) +
        geom_point() +
        geom_smooth(se = FALSE)+
        ggtitle("AvgAirDose against Distance to Daichi Plant")

1.19 Nihommatsu 2011,Counts of Measuring Locations per Air Dose Rate
ggplot(data = fuk2013_q) +
        geom_bar(mapping = aes(x = daichi_distance, fill = dose_quants), width = 1)+
        ggtitle("AvgAirDose Measured Counts against Daichi Distance")

2.10 Air Dose Rate Reduction due to Distance to Fukushima Daichi

\[AvgAirDoseRate = \beta_{0} + \beta_{1} daichi distance + \epsilon_{i}\]

fit1 <- lm(AvgAirDoseRate~daichi_distance, data = fuk2011)
summary(fit1)
## 
## Call:
## lm(formula = AvgAirDoseRate ~ daichi_distance, data = fuk2011)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.226 -0.901 -0.321  0.165 34.291 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      2.7838704  0.0267272  104.16   <2e-16 ***
## daichi_distance -0.0288916  0.0004048  -71.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.278 on 45271 degrees of freedom
## Multiple R-squared:  0.1012, Adjusted R-squared:  0.1011 
## F-statistic:  5095 on 1 and 45271 DF,  p-value: < 2.2e-16
#Confidence interval
confint(fit1)
##                       2.5 %      97.5 %
## (Intercept)      2.73148462  2.83625613
## daichi_distance -0.02968495 -0.02809831

2.20 Air Dose Rate Reduction due to Population Density and Land Use

\[AvgAirDoseRate = \beta_{0} + \beta_{1} popn density + \beta_{2} human activities + \epsilon_{i}\]

Challenge

The Air Dose Rate measurements were collected from three sources:

  1. Cars and Buses with Sensors: Measure Air of 100m&#0178 from the road
  2. Un manned Auto Helicopter: Measure from 300m above ground
  3. Fixed Sensors: Geographically sparsed

What constitutes Land Use?:

  1. Farming: Short life span crops lead to frequent land tilling
  2. Cleaned Places: Parks, GolfCourses, Gardens are cleaned too and could reduce radiations

Population density and land usage is measured on a 500m² basis

fuk_pop <- read.csv("fuk.csv")
# Convert mesh gridecode to lat/long
mesh_latlong <- function (mesh, loc = "C")
{
    mesh <- as.character(mesh)

    if (length(grep("^[0-9]{4}", mesh)) == 1) { 
        mesh12 <- as.numeric(substring(mesh, 1, 2))
        mesh34 <- as.numeric(substring(mesh, 3, 4))
        lat_width  <- 2 / 3;
        long_width <- 1;
    }
    else {
        return(NULL)
    }

    if (length(grep("^[0-9]{6}", mesh)) == 1) { 
        mesh5 <- as.numeric(substring(mesh, 5, 5))
        mesh6 <- as.numeric(substring(mesh, 6, 6))
        lat_width  <- lat_width / 8;
        long_width <- long_width / 8;
    }

    if (length(grep("^[0-9]{8}", mesh)) == 1) { 
        mesh7 <- as.numeric(substring(mesh, 7, 7))
        mesh8 <- as.numeric(substring(mesh, 8, 8))
        lat_width  <- lat_width / 10;
        long_width <- long_width / 10;
    }

    lat  <- mesh12 * 2 / 3;          
    long <- mesh34 + 100;

    if (exists("mesh5") && exists("mesh6")) {  
        lat  <- lat  + (mesh5 * 2 / 3) / 8;
        long <- long +  mesh6 / 8;
    }
    if (exists("mesh7") && exists("mesh8")) {  
        lat  <- lat  + (mesh7 * 2 / 3) / 8 / 10;
        long <- long +  mesh8 / 8 / 10;
    }

    if (loc == "C") {   
        lat  <-  lat  + lat_width  / 2;
        long <-  long + long_width / 2;
    }
    if (length(grep("N", loc)) == 1) {  
        lat  <- lat  + lat_width;
    }
    if (length(grep("E", loc) == 1)) {  
        long <- long +long_width;
    }

    lat  <- sprintf("%.8f", lat); 
    long <- sprintf("%.8f", long);

    x <- list(as.numeric(lat), as.numeric(long))
    names(x) <- c("lat", "long")

    return (x)
}
# j <- mesh_latlong(mesh = 564000463 , loc = "C")
# class(j);print(j);length(j)
grides <- fuk_pop[,1]
head(grides);class(grides);length(grides)
## [1] 564000194 564000364 564000382 564000393 564000462 564000463
## [1] "integer"
## [1] 10831
#create the lat/long
mylist <- list()
for (i in 1:length(grides) ){
        lis <- mesh_latlong(mesh = grides[i], loc = "C")
        mylist[[i]] <- lis
}

res <- as.data.frame(mylist)


#test
# df <- data.frame(matrix(unlist(res), nrow=10831, byrow=T))
# df <- ldply (res, data.frame)

library(data.table)
df <- as.data.frame(rbindlist(mylist, fill=TRUE))
df[,"gridcode"] <- grides 

# merge gridcode and lat/lon datasets of Fuk population
fuk_ll <- merge(fuk_pop, df, by.x = "gridcode", by.y = "gridcode", all = TRUE)
write.csv(fuk_pop, file="fuk_pop.csv",row.names = FALSE)
#load new dataset
fuk_pop <- read.csv("fuk_pop.csv")
View(fuk_pop)
no_unique <- length(unique(fuk_pop$gridcode))
no_unique
## [1] 10831

trial plots

pop_plot <- leaflet() %>%
        addTiles()%>%
        addPolylines(data = fuk_ll,lng = ~long, lat = ~lat)
pop_plot        
#remove duplicate lat/lon
uniq <- unique(x=df[,-3])

uniq_plot <- leaflet() %>%
        addTiles()%>%
        addCircles(data = uniq,lng = ~long, lat = ~lat)
uniq_plot
# Unique, lat or lon only
data <- subset(fuk_ll, !duplicated(fuk_ll[,7]))
dim(data)
## [1] 80  7
uniq_plot <- leaflet() %>%
        addTiles()%>%
        addCircles(data = data,lng = ~long, lat = ~lat)
uniq_plot

Understand population distribution and airdose

Population link:http://www.e-stat.go.jp/SG1/estat/toukeiChiri.do?method=init Fuk Population 5640: 1degree, 40min,80km 00: 5min,7’30min, 10km 19: 30s, 45s, 1km 4: 15s, 22.5s, 500m [1:SW,2:SE,3:NW,4:NE]

#remove slashes from 2013 AirDose grids
fuk2013_q$gridcode <- gsub("_","",fuk2013_q$gridcode)
#change fuk_ll$gridcode to character from integer
fuk_ll$gridcode <- as.character(fuk_ll$gridcode)
# Merge fuk_ll and fuk2013_q

fuk_air_pop <- merge(fuk2013dose,fuk_pop, by.y = "gridcode", all.x )

2.30 Using other machine learning algorithms

The Fukushima Nuclear Radiations is multi dimensional;

  • Three major Isotopes each with;
  • Time series due to half life
  • Magnitude of radiation (Becquerel)
  • Absorbed dose (Sievert (Sv))
  • The above dimensionalities coupled with distance,population density and land use, create a data set that can be run on extensive machine learning algorithms like Support Vector Machines, Random Forest, Recurrent Neural Networks and more.

    end

    2.40 LAND USE

    http://nlftp.mlit.go.jp/ksj-e/old/cgi-bin/_download_view.cgi Forest and Public Lands mesh Land Classification Mesh

    Agricultural Areas (surface):A12-60A-48-01.1.zip 10.8MB Agricultural Areas (surface):A12-60A-48-01.0a.zip 10.82MB

    Airdose change to 500m2